function [] = computeElasticities_providerExits(runId)

%%% Add paths
addpath('../model_funcs');
addpath('../general_funcs');

%%% Get input / output paths
config = getConfig_providerExits(runId);

if config.modelId ~= 3
	error('This code is only for modelId == 3 (continuous time models).');
end

outputFolder = createFolderIfNotExist(sprintf('%s/elasticities', config.runFolder));
outputMatFile = sprintf('%s/elasticities.mat', outputFolder);
outputCSVFile1  = sprintf('%s/elasticities.csv', outputFolder);
outputCSVFile2  = sprintf('%s/elasticitiesStrings.csv', outputFolder);

%%% Load data
load(sprintf('%s/runData.mat', config.runFolder), 'data');
%% Get paramNames
paramNames = get_param_names_providerExits(data);

%%%  Load parameter estimates from the 3 submodels
load(sprintf('%s/results.mat', config.runFolder), 'params_hat', 'params_ses');
params = params_hat;
%%% Split params_parts + save everything in object paramsParts
paramsParts  = split_combine_parts(data.dims,  1, params);

%%% Read dimensions
NumSpells = data.dims.dims1{2};
K = data.dims.dims1{3};
NumParams = data.dims.NumParams;
NumParts = length(data.Xparts);
	
myfunc1 = @(x) x; % Identity function
myfunc2 = @(x) 1; % one
myfunc3 = @(x) 1-exp(-x); % log(1+x) function

elastNames = {'wrt_X', 'wrt_eX', 'wrt_eX_minus1', 'wrt_Area_TIMES_eX_minus1'};


%%%% IDEA: for each spell, compute average elasticities across event types (adoptions in canton k), then take a weighted average across
% spells, where the weight is the duration of the spell.

elasticities1 = zeros(NumParams, NumSpells); % Need real need to store all spell results...
elasticities2 = zeros(NumParams, NumSpells); % Need real need to store all spell results...
elasticities3 = zeros(NumParams, NumSpells); % Need real need to store all spell results...

% Assign things that won't depend on spell
data_uu.surfaces                 = data.surfaces;
data_uu.dims.dimType             = data.dims.dimType;
data_uu.dims.NumParts            = data.dims.NumParts;
data_uu.dims.NumParams           = data.dims.NumParams;
data_uu.dims.Xpart_2_NumX        = data.dims.Xpart_2_NumX;
data_uu.dims.Xpart_2_NumX_FEs    = data.dims.Xpart_2_NumX_FEs;
data_uu.dims.Xpart_2_Num_FE_vals = data.dims.Xpart_2_Num_FE_vals;
data_uu.dims.NumFEvals2Keep      = data.dims.NumFEvals2Keep;
data_uu.dims.NumObs   = K;
data_uu.dims.dims1{1} = K;
data_uu.dims.dims1{2} = 1;
data_uu.dims.dims1{3} = K;

for pp = 1:NumParts
	data_uu.Xparts{pp}.NumX_FE_vals = data.Xparts{pp}.NumX_FE_vals;
	data_uu.Xparts{pp}.Xnames       = data.Xparts{pp}.Xnames;
	data_uu.Xparts{pp}.X_FE_labels  = data.Xparts{pp}.X_FE_labels;
end
data_uu.Xparts{3} = data.Xparts{3};

X1 = data.Xparts{1}.X; % (NumSpells*K) x NumX1
X1 = reshape(X1, [NumSpells K size(X1,2)]); % NumSpells x K x NumX1
X1 = permute(X1, [2 3 1]); % K x NumX1 x NumSpells

X_FEs1 = data.Xparts{1}.X_FEs; % (NumSpells*K) x NumX_FEs1
X_FEs1 = reshape(X_FEs1, [NumSpells K size(X_FEs1,2)]); % NumSpells x K x NumX_FEs1
X_FEs1 = permute(X_FEs1, [2 3 1]); % K x NumX_FEs1 x NumSpells

denom = data.surfaces; % K x 1

tic
for uu = 1:NumSpells
	data_uu.Xparts{1}.X     = X1(:,:,uu);     % K x NumX1
	data_uu.Xparts{1}.X_FEs = X_FEs1(:,:,uu); % K x NumX_FEs1
	
	data_uu.Xparts{2}.X     = data.Xparts{2}.X(uu,:); % 1 x NumX
	data_uu.Xparts{2}.X_FEs = data.Xparts{2}.X_FEs(uu,:); % 1 x NumX_FEs

	elasts1_uu = compute_mean_Xfunc_times_beta(data_uu.dims, data_uu.Xparts, myfunc1, paramsParts)';
	elasts2_uu = compute_mean_Xfunc_times_beta(data_uu.dims, data_uu.Xparts, myfunc2, paramsParts)';
	elasts3_uu = compute_mean_Xfunc_times_beta(data_uu.dims, data_uu.Xparts, myfunc3, paramsParts)';

	% Compute elasts4 wrt Area_TIMES_eX_minus1 (for variables X = log([1+Xtilda]/Area) )
	Xfunc_mean_beta_parts = cell(NumParts,1);
	for pp = 1:NumParts
		Xpart_pp = data_uu.Xparts{pp};
		if pp == 1 % T-K part
			NumXp = size(Xpart_pp.X, 2);
			X_pp = reshape(Xpart_pp.X, [1 K NumXp]);
			Xfunc_mean_beta_parts{pp}.X_Y = paramsParts{pp}.beta .* mean(reshape(1 - exp(-X_pp)./denom', [K NumXp]), 1)'; % NumXp x 1
		end
		if pp == 2 % T-part
			Xfunc_mean_beta_parts{pp}.X_Y = NaN * ones(data_uu.dims.Xpart_2_NumX(pp),1); % NumX x 1
		end
		if pp == 3 % K-part
			Xfunc_mean_beta_parts{pp}.X_Y = paramsParts{pp}.beta .* mean(1 - exp(-Xpart_pp.X)./denom, 1)'; % NumX x 1 
		end
		NumX_FEs = size(Xpart_pp.X_FEs,2);
		Xfunc_mean_beta_parts{pp}.X_FEs_Y = cell(NumX_FEs,1);
		for ff = 1:NumX_FEs
			NumFE_vals_ff = Xpart_pp.NumX_FE_vals(ff);
			Xfunc_mean_beta_parts{pp}.X_FEs_Y{ff} = NaN * ones(NumFE_vals_ff,1);
		end
	end
	elasts4_uu = split_combine_parts(data_uu.dims, 3, Xfunc_mean_beta_parts)'; % 1 x NumParams
	
	% Store everything
	elasticities1(:,uu) = elasts1_uu;
	elasticities2(:,uu) = elasts2_uu;
	elasticities3(:,uu) = elasts3_uu;
	elasticities4(:,uu) = elasts4_uu;

	displayRemainingTime(uu, NumSpells, 1000);
end

spellWeights = data.spellDurations/sum(data.spellDurations); % NumSpells x 1
elasts1 = elasticities1*spellWeights;
elasts2 = elasticities2*spellWeights;
elasts3 = elasticities3*spellWeights;
elasts4 = elasticities4*spellWeights;

allElasticities = [elasts1 elasts2 elasts3 elasts4]; % NumParams x 4

% Compute standard errors on elasticities
allElasticities_ses = allElasticities .* params_ses./params;


%%%% Get the relevant elasticity for each variable, depending on any transformation applied
paramNames = strrep(paramNames, 'Beta ', '');
paramNames = strrep(paramNames, ' on ProviderExits', '');
rowNamesMapping = getRowNamesMapping();
[flag,idxes] = ismember(paramNames, rowNamesMapping(:,1));
idxes                = idxes(flag);
allElasticities      = allElasticities(flag,:);
allElasticities_ses  = allElasticities_ses(flag,:);
paramNames           = rowNamesMapping(idxes,2);

myVariablesAndElastType = {
	'Log density of providers at destination at k1'               , 4;
	'Log density of consumers living in destination at k1'        , 4;
	'Log density of consumers who transacted in destination at k1', 4;
	'Log density of population at k1'                             , 2;
	'Log density of incoming commuters at k1'                     , 4;
	'Log density of tourist beds at k1'                           , 4;
	'Costs related to car providership (in 1,000 euros) at k1'    , 1;
	'Income (in 100,000 euros) at k1'                             , 1;
	'Unemployment rate at k1'                                     , 1;
	'Frac. of people without high school diploma at k1'           , 1;
	'Frac. of people with higher educ. diploma at k1'             , 1;
	'Fraction of votes for ecological candidate at k1'            , 1; 
	'Google Trend of competitor at k1'                            , 1; 
	'Precipitation at destination (in decimeters) at k1'          , 1;
	'Temperature at destination (in celsius) at k1'               , 1;
};
[flag2,idxes] = ismember(myVariablesAndElastType(:,1), paramNames); assert(all(flag2));
paramNames = paramNames(idxes);
allElasticities     = allElasticities(idxes,:);
allElasticities_ses = allElasticities_ses(idxes,:);
NumElasts = length(myVariablesAndElastType);
myelasts = zeros(NumElasts,1);
myelasts_ses = zeros(NumElasts,1);
for ee = 1:NumElasts
	myelasts(ee)     = allElasticities(ee,myVariablesAndElastType{ee,2});
	myelasts_ses(ee) = allElasticities_ses(ee,myVariablesAndElastType{ee,2});
end

% Rename those paramNames with log-density (after the fact)
paramNames = strrep(paramNames, 'Log density', 'Density');

M1 = table(myelasts, myelasts_ses);
M1.Properties.RowNames = paramNames;
writetable(M1, outputCSVFile1, 'WriteRowNames', true);

M2 = cell(NumElasts, 1);
for ee = 1:NumElasts
	M2{ee} = sprintf('%.2f (%.2f)', myelasts(ee), myelasts_ses(ee));
end
M2 = cell2table(M2);
M2.Properties.VariableNames = {'Elasticity'};
M2.Properties.RowNames = paramNames;
disp(M2);
writetable(M2, outputCSVFile2, 'WriteRowNames', true);

end
